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1.  Introduction 


An  accurate  understanding  of  fundamental  flow  physics  and  accurate  prediction  of 
aerodynamic  coefficients  is  critical  to  the  development  of  precision  munitions 
(Fresconi  et  al.  2011;  Coyle  and  Silton  2015;  Silton  and  Fresconi  2015).  The  flow 
fields  associated  with  precision  munitions  can  be  very  complex,  involving  shock¬ 
boundary  layer  interaction,  unsteady  wakes,  and  body,  canard  and  tail-fin 
interactions  (Sahu  1990;  DeSpirito  et  al.  2003;  Silton  2005;  DeBonis  et  al.  2012; 
Bhagwandin  2014,  2015;  Scheuermann  et  al  2015;  Silton  and  Coyle  2015;  Silton 
and  Fresconi  2015).  Simple  empirical  and  analytical  techniques  such  as  PRODAS 
(Arrow  Tech  Associates  2015),  Missile  Datcom  (Rosema  et  al.  2011),  MISL3 
(Lesieutre  et  al.  2002),  and  AP09  (Moore  and  Moore  2008)  are  not  always  able  to 
accurately  determine  the  complex  3-D  flow  interactions  and  associated  nonlinear 
flow  physics.  Computational  Fluid  Dynamics  (CFD)  offers  an  alternative  to 
compute  these  nonlinear  interacting  flow  fields,  providing  detailed  understanding 
of  the  associated  nonlinear  aerodynamic  processes  and  predicting  the  associated 
aerodynamic  coefficients  that  are  required  to  determine  if  performance 
requirements  can  be  met.  CFD  has  continued  to  emerge  as  a  critical  technology  for 
the  aerodynamic  design  and  assessment  of  flight  vehicles.  Bringing  CFD  into  the 
earlier  stages  of  design  and  development  is  necessary  for  it  to  be  more  effective;  to 
bring  CFD  into  the  earlier  stages,  though,  the  results  must  be  obtained  as  quickly 
and  accurately  as  possible.  However,  even  with  the  advancement  of  high- 
performance  computing,  CFD  computing  time  is  still  much  greater  than  the  time 
required  to  obtain  solutions  using  engineering-level  codes.  This  means  that  CFD  is 
not  implemented  until  further  along  in  the  process,  often  after  a  problem  is 
discovered. 

Steady-state  and/or  time-accurate  CFD  methods  can  be  used  to  obtain  a  complete 
set  of  aerodynamic  coefficients  (static  coefficients  and  dynamic  derivatives)  in  all 
flight  regimes,  from  subsonic  to  supersonic  speeds.  In  general,  most  of  the 
aerodynamic  coefficients  of  a  projectile  or  missile  in  the  supersonic  regime  can  be 
generated  efficiently  using  steady-state  methods;  only  pitch  damping  may  need  to 
be  calculated  in  a  time-accurate  manner.  Time-accurate  methods  are  sometimes 
required  to  accurately  predict  base  flows  and  dynamic  derivatives  of  these 
munitions,  especially  in  the  subsonic  and  transonic  flight  regimes  (Oktay  and  Akay 
2002;  Park  and  Kwon  2003;  DeSpirito  and  Heavey  2004;  Murman  2007;  DeSpirito 
and  Plostins  2007;  DeSpirito  2008;  DeSpirito  et  al.  2009;  Silton  201 1;  Da  Ronch  et 
al.  2012;  Bhagwandin  2012;  Klatt  et  al.  2012;  Bhagwandin  and  Sahu  2014).  The 
static  aerodynamic  coefficients  can  be  obtained  from  a  steady-state  method.  A 
number  of  researchers  (Weinacht  and  Sturek  1990;  DeSpirito  et  al.  2009; 
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Bhagwandin  and  Sahu  2014)  found  the  pitch-damping  coefficient  derivative  could 
also  be  determined  using  a  steady-state  coning  methodology.  More  recent  advances 
in  computational  modeling  have  also  led  to  coupling  of  CFD  codes  to  rigid-body 
dynamics  (RBD)  codes  for  the  time-accurate  simulation  of  free-flight  motion, 
which  allows  determination  of  both  the  static  aerodynamic  coefficients  and 
dynamic  derivatives  using  the  same  numerical  simulation  (Sahu  2004;  Costello  et 
al.  2007;  Sahu  2008;  Sahu  2009;  Stahl  et  al.  2009;  Montalvo  and  Costello  2010; 
Sahu  2011;  Camargo  et  al.  2012;  Sahu  and  Fresconi  2015;  Wang  et  al.  2015). 

When  designing  a  new  precision  munition,  especially  in  the  early  stages,  it  is  the 
static  aerodynamic  coefficients  that  are  of  greatest  importance  to  ensure  the 
munition  can  meet  its  range  and  maneuverability  requirements.  This  requires  the 
development  of  an  aerodynamic  database  that  likely  includes  a  range  of  angles  of 
attack  and  either  roll  angles  or  side-slip  angles  at  a  given  (or  multiple)  Mach 
number(s).  While  the  new  coupled  CFD-RBD  method  would  require  fewer  total 
simulations,  it  is  quite  time  intensive  as  it  must  be  solved  using  the  time-dependent 
Reynolds-Averaged  Navier-Stokes  (RANS)  equations  with  very  small  time  steps 
to  accurately  capture  the  flow  physics.  The  traditional  steady-state  solution  uses  the 
same  governing  RANS  equations,  but  allows  for  larger  time  steps  to  more  quickly 
converge  the  solution.  Flowever,  a  large  number  of  these  steady-state  solutions 
would  be  required  to  generate  the  aerodynamic  database  when  nonlinear 
aerodynamic  behavior  is  present  (canard  stall,  vortex  interactions,  etc.)  in  the  range 
of  angle  of  attack  of  interest  causing  this  method  to  be  time  intensive  as  well. 

Recently,  Sahu  and  Heavey  (2010)  demonstrated  an  alternative  method  for  rapid 
determination  of  static  aerodynamic  coefficients  using  the  time-accurate  RANS 
equations;  a  quasi-steady  sweep  procedure  that  allows  determination  of  the 
coefficients  for  an  entire  range  of  angles  of  attack.  Although  their  research  was  on 
a  complex  nonaxisymmetric  missile  configuration,  it  was  limited  to  the  transonic 
and  supersonic  flight  regimes.  The  work  presented  in  this  paper  extends  their  work 
to  the  subsonic  regime.  Numerical  computations  using  the  quasi-steady  sweep 
procedure  are  performed  at  Mach  0.65  on  the  high-maneuverability  airframe 
(HMA)  being  developed  at  the  US  Army  Research  Laboratory  (ARL)  as  a 
demonstrator  platform  for  the  precision-munitions  research  area  (Fresconi  et  al. 
2011;  Silton  and  Fresconi  2013;  Fresconi  et  al.  2014;  Silton  et  al.  2014). 
Steady-state  simulations  were  also  performed  for  comparison  with  the  sweep 
results.  Computational  resource  requirements  are  included  to  show  the  efficiency 
of  the  quasi-steady  sweep  procedure. 
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2.  Solution  Technique 


The  commercially  available  code  CFD 1 1  vl2.1.1  (Metacomp  Technologies  201 1), 
by  Metacomp  Technologies,  Inc.,  is  used  for  the  CFD  simulations.  CFD 1 1  is  a 
finite-volume,  unstructured  solver  capable  of  computing  CFD  solutions  for  a  wide 
range  of  aerospace  applications.  In  this  study,  CFD++  is  used  in  both  its  steady-state 
and  time-accurate  modes  to  numerically  solve  the  3-D,  compressible  RANS 
equations  to  compute  the  flow  solution. 

The  Goldberg  3 -equation  k-s-Rt  turbulence  model  (Goldberg  et.  al.  1998)  was 
implemented  in  this  study  based  on  the  findings  by  Silton  and  Fresconi  (20 13).  This 
3-equation  model  solves  the  transport  equations  for  undamped  eddy  viscosity,  Rt, 
in  addition  to  turbulent  kinetic  energy  ( k )  and  its  dissipation  rate  (s).  This  solution 
methodology  accounts  for  nonequilibrium  conditions  and  avoids  free-stream- 
turbulence  decay  under  shear-free  flow  conditions.  For  the  current  problem, 
initialization  of  the  turbulence  transport  was  completed  by  setting  the  turbulence 
intensity  to  2%  and  the  turbulent-to-molecular  viscosity  ratio  to  50,  as  the  length 
scale  was  not  known. 

2.1  Steady-State  Simulation 

The  flow  solution  was  advanced  toward  steady-state  convergence  using  a 
point-implicit  time-integration  scheme  with  local  time-stepping,  defined  by  the 
Courant-Friedrichs-Lewy  (CFL)  number.  A  linear  ramping  schedule  was  used  to 
gradually  increase  the  CFL  number  over  the  first  1 00  iterations  after  which  the  CFL 
remained  constant  until  convergence  was  reached.  As  only  the  subsonic  flight 
regime  (0.3  <  M<  0.8)  is  considered  in  this  paper,  the  suggested  CFL  ramping  is 
from  1  to  100.  Depending  on  the  angle  of  attack  being  investigated,  the  maximum 
CFL  number  for  the  ramping  was  limited  to  25  or  50.  The  multigrid  W-cycle 
method  with  a  maximum  of  4  cycles  and  20  coarse  grid  levels  was  used  to 
accelerate  convergence.  Implicit  temporal  smoothing  was  applied  for  increased 
stability,  which  is  especially  useful  where  strong  transients  arise.  The  spatial 
discretization  function  was  a  second-order,  upwind  scheme  using  a 
Harten-Lax-van  Leer-Contact  (commonly,  HLLC)  Riemann  solver  and 
Metacomp’s  multidimensional  Total-Variation-Diminishing  (commonly,  TVD) 
flux  limiter  (Metacomp  Technologies  2011).  Double-precision  format  was  used  for 
all  computations. 

Solution  convergence  (typically  within  3,000  iterations)  was  determined  by  1)  a 
several-orders-of-magnitude  reduction  of  the  cell -averaged  residuals  of  the  RANS 
equations  until  the  residuals  plateaued  and  2)  the  total,  tail-fin,  and  canard  forces 
and  moments  reaching  asymptote  and  unchanging  to  within  0.1%. 
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2.2  Quasi-Steady  Sweep  Simulation 


Using  the  equivalent  of  a  series  of  steady-state  simulations,  the  flow  solution  for 
the  quasi-steady  sweep  is  generated  using  the  CFD  solver’s  time-accurate  mode  to 
obtain  the  static  aerodynamic  coefficients  across  a  range  of  angles  of  attack  in  a 
single  numerical  simulation.  The  time-accurate  mode  utilizes  the  implicit  solver 
with  dual  time-stepping.  Dual  time-stepping  employees  2  time-steps;  an  “outer”  or 
global  time-step  that  corresponds  to  the  time  discretization  of  the  physical  time 
variation  and  an  “inner”  or  local  time  variation  that  helps  to  satisfy  the  physical 
transient  equations. 

The  “outer”  time-step  is  chosen  to  provide  the  desired  time  accuracy  and  is  applied 
to  every  cell.  Typically,  this  is  set  to  a  small  value  to  ensure  that  the  flow  features 
are  sufficiently  resolved  in  time.  In  the  present,  sweep  simulation,  the  time- 
accurate,  dual  time-step  method  is  used  to  determine  a  “quasi-steady”  solution  at 
each  time-step,  which  corresponds  to  an  increment  of  rotation,  in  an  efficient 
manner.  Therefore,  a  large  global  time-step  of  1 .0  s  is  employed.  This  large  time- 
step  drives  the  rate  of  change  of  the  flow  variables  with  respect  to  angle  of  attack 
to  zero  as  quickly  as  possible.  Thus,  the  solutions  are  not  transient,  as  is  typical  of 
a  time-accurate  simulation,  but  rather  multiple  static  solutions. 

The  “inner”  time-step  is  determined  internal  to  the  solver  via  the  CFL  number  that 
is  equal  to  the  final  maximum  CFL  number  from  the  steady-state  solution.  The 
inner  time-step  (or  iteration)  is  a  local,  nonphysical  time-step  used  to  converge  the 
time-accurate  RANS  equations  at  each  physical  time-step.  The  time-step  for  the 
inner  iteration  is  allowed  to  vary  spatially.  For  the  inner  iteration,  the  multigrid 
W-cycle  method  and  implicit  temporal  smoothing  are  applied  to  accelerate 
convergence  of  the  inner  iterations  and  therefore  the  outer  global  time-step.  If  the 
inner  iterations  are  not  sufficiently  converged,  the  solution  is  not  converged  and 
there  is  a  loss  in  accuracy.  The  number  of  inner  iterations  is  chosen  such  that 
approximately  a  2-order-of-magnitude  reduction  in  the  cell-averaged  “inner” 
residuals  was  achieved.  Typically,  100  inner  iterations  are  found  to  be  adequate  in 
the  present  simulations  using  the  quasi-steady  sweep  procedure.  These 
requirements  are  significantly  larger  than  for  a  “standard”  time-accurate  solution 
where  only  5-20  inner  iterations  are  required  for  a  1 -2-order-of-magnitude  drop  in 
the  cell-averaged  inner  residuals  and  a  converged  solution. 

The  steady-state  solution  with  the  mesh  rotated  to  the  lowest  desired  angle  of  attack 
serves  as  the  initial  condition  for  the  initial  sweep  solution.  For  each  step  of  the 
sweep  solution,  the  computational  grid  is  rotated  by  a  small  increment  and  a  new 
solution  is  computed.  In  CFD++,  the  rotation  takes  place  at  the  end  of  the  time-step; 
the  first  time-step  is  solved  at  the  starting  angle  of  attack.  This  process  continues 
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across  the  desired  range  of  angles  of  attack.  Figure  1  shows  the  location  of  the  body 
at  various  stages  of  the  sweep.  In  this  example,  the  rotation  point  is  at  the 
projectile’s  center  of  gravity.  The  rotation  angle  must  be  chosen  small  enough  to 
capture  any  nonlinear  aerodynamic  effects  that  may  be  present.  In  the  present  sweep 
simulations,  a  rotation-angle  increment  of  0.25°  was  able  to  suitably  capture  the 
nonlinear  aerodynamic  effects. 


Fig.  1  Location  of  the  projectile  at  various  stages  of  the  sweep  procedure 

Every  numerical  solver  has  its  nuances  in  how  the  sweep  procedure  is  implemented. 
Details  of  implementation  of  the  quasi-steady  sweep  simulations  in  CFD 1 1  can  be 
found  in  the  Appendix. 

3.  Geometry 


3.1  Model  Geometry 

The  geometry  utilized  for  this  study  was  the  HMA  designed  as  the  platform  to 
demonstrate  advanced  guided  munitions  using  low-cost  technologies.  Both  the 
body-fin  configuration  (BF)  and  the  body-fin-canard  configuration  were 
considered  (Fig.  2). 
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Fig.  2  Subsonic  geometries:  BF  (top)  and  body-fin-canard  configuration  (bottom) 

The  geometric  model  was  created  in  SolidWorks  (Dassault  Systemes  2012)  and  has 
a  body  reference  diameter  of  83  mm  (1  caliber).  This  was  a  simplified  geometry 
with  the  gaps  and  slots  in  the  body  needed  to  accommodate  the  tail  fins  suppressed. 
The  body  geometry  consisted  of  a  hemispherical  nose  cap;  a  nearly  cylindrical  body 
320.38  mm  (3.86  cal.)  long;  and  a  66.4-mm  (0.8 -cal) -long,  7°  boattail.  The  overall 
projectile  length  was  427.23  mm  (5.15  cal.).  Ten  flat-plate  fins  having  a  span  of 
217.22  mm  and  an  average  chord  of  20.5  mm,  with  the  trailing  edge  located 
4.48  mm  from  the  base  of  the  projectile,  were  included  for  stability.  The  tail  fins 
were  shaped  such  that  they  can  store  flush  in  the  body  for  gun  launch  while 
maximizing  area  for  stability.  A  7-mm-high,  3.8-mm-deep  tab  at  the  tip  of  each  tail 
fin  was  bent  at  a  15°  angle  such  that  the  projectile  would  roll  clockwise  (looking 
from  rear  of  airframe). 

For  the  body-fin-canard  configuration,  4  NACA0015  airfoils,  rotated  45°  off 
plane,  were  added  to  the  model.  Each  of  the  4  canards  had  a  chord  of  18.86  mm 
and  a  semispan  of  96.36  mm  with  the  quarter  chord  located  144.7  mm  aft  of  the 
nose  of  the  projectile.  Canard  deflection  angles  of  0°  and  4°  pitch  up  about  the 
quarter  chord  were  considered.  For  simplicity,  the  model  had  the  canards  directly 
attached  to  the  body;  the  body  slots  were  suppressed. 


3.2  Numerical  Grids 
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The  grids  used  in  the  numerical  simulation  were  created  with  MIME  v4.1  by 
Metacomp  Technologies  (Metacomp  Technologies  2010).  MIME  is  an 
unstructured  mesh  generator  that  allows  triangular  or  quadrilateral  dominant  cells 
for  the  surface  mesh.  Once  an  adequate  surface  mesh  is  generated,  prism  layers  can 
be  specified  and  created  when  the  volume  mesh  (tetrahedral  dominant)  is  generated. 

Figure  3  shows  the  extent  of  the  outer  boundary.  It  extended  approximately  3  5  body 
lengths  from  the  projectile  in  all  directions.  Cylindrical  density  boxes  were  placed 
in  the  wake  of  the  projectile  and  tail  fins  to  ensure  proper  resolution  in  this  area. 
An  additional  density  box  (No.  3)  was  placed  between  the  canards  and  the  tail  fins, 
when  the  canards  were  present,  to  ensure  the  canard-tip  vortices  would  not  dissipate 
significantly  due  to  mesh  resolution  prior  to  reaching  the  tail  fins.  The  surface  mesh 
and  boundary-layer  growth  was  restricted  to  a  ratio  of  1 .2  or  less.  The  growth  ratio 
was  relaxed  to  2.0  when  creating  the  remainder  of  the  volume  mesh.  The  prism- 
layer  spacing  was  chosen  such  that  wall-function  spacing  (nondimensional  wall 
distance,  30  <  y  +  <  60)  was  present  on  the  body,  while  a  solve-to-wall  mesh 
y+  w  1  was  used  on  and  around  the  tail  fins  and  canards.  To  achieve  this,  a  first 
cell  spacing  for  the  prism  layer  was  specified  at  3.5  x  10-3  mm  on  the  tail  fins  and 
canards,  0.25  mm  on  the  body,  and  0.2  mm  on  the  nose.  This  produced  y+  *  1  on 
the  tail  fins  and  canards  as  well  as  on  the  body  immediately  surrounding  them.  On 
the  nose  of  the  airframe  and  remainder  of  the  body,  the  y  +  varied  between  25  and 
45.  MIME  allowed  for  a  smooth  transition  between  the  specified  spacings.  Mesh- 
refinement  studies,  which  showed  that  this  mesh  achieved  grid  independence,  have 
been  conducted  for  this  body-fin  geometry  (Silton  and  Fresconi  2013)  and  are  not 
included  in  this  paper.  Mesh-refinement  studies  for  the  body-fin-canard  geometry 
have  not  yet  been  completed.  The  computational  domain  is  fully  3-D  with  no 
symmetry.  The  body-fin  mesh  consisted  of  approximately  26  million  cells,  while 
the  body-fin-canard  mesh  consisted  of  approximately  45  million  cells. 
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Fig.  3  Extent  of  outer  mesh  boundary  with  close-up  of  projectile  showing  locations  of 
density  boxes 

3.3  Boundary  Conditions 

The  entire  far-field  boundary  was  set  as  “characteristics-based”  inflow/outflow. 
This  boundary  condition  takes  the  specified  free-stream  conditions  and  solves  a 
Riemann  problem  at  the  boundary  using  the  supplied  data  as  a  virtual  state  outside 
the  domain.  The  walls  of  the  projectile  were  specified  as  an  adiabatic,  no-slip, 
viscous  boundary.  Wall  functions  were  used  to  compute  the  boundary  layer  on  the 
projectile’ s  surface  except  where  there  was  adequate  prism-layer  resolution  to  solve 
to  the  wall,  as  was  the  case  on  the  tail  fins  and  canards. 

4.  Results  and  Discussion 


Numerical  simulations  were  performed  to  assess  the  ability  and  efficiency  of  the 
quasi-steady  sweep  procedure  to  predict  the  flow  field  and  aerodynamic 
coefficients  for  a  complex  gun-launched-projectile  configuration  in  the  subsonic 
flight  regime.  All  computations  were  completed  using  a  free-stream  temperature 
and  pressure  of  288.15  K  and  101325  Pa,  respectively.  The  study  was  conducted  at 
Mach  0.65  (220  m/s).  Angles  of  attack  between  -14°  and  14°,  inclusive,  were 
considered.  The  flow  domain  was  initialized  using  free-stream  conditions 
everywhere. 
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4.1  Body-Fin  Configuration 


Steady-state  simulations  were  completed  for  the  BF  at  a  =  -14°,  0°,  2°,  3°,  5°,  8°, 
10°,  and  14°.  The  quasi-steady  sweep  procedure  was  performed  using  the  steady- 
state  solution  at  a  =  —14°  as  the  initial  condition.  The  critical  parameters  were 
varied  to  determine  the  requirements  for  an  accurate  solution.  The  critical 
parameters  varied  were  CFL  =  10  and  50;  angle-of-attack  increments  (da)  =  0.5°, 
0.25°,  and  0.1  °;  and  number  of  inner  iterations  (N)  =  50,  100,  and  200. 

4.1.1  CFL  =10 

Initially,  parameters  were  set  to  CFL  =10  (small  CFL  is  typically  warranted  when 
determining  time-accurate  solutions)  with  50  inner  iterations  (Sahu  and  Heavey 
2010)  and  da  =  0.25°  and  0.1°  were  investigated.  Figure  4  shows  the  convergence 
history  of  the  solution  residual  for  the  inner  iterations  as  a  function  of  the  da  used. 
Both  increments  produce  approximately  1  -order-of-magnitude  drop  in  residuals, 
although  they  do  not  level  off  as  would  be  expected  of  a  completely  converged 
solution. 
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Residual  Plot  (Normalized) 


Residual  Plot  (Normalized) 


Fig.  4  Time  history  of  solution  residual  as  a  function  of  da  for  CFL  =  10,  N  =  50; 
da  =  0.25°  (top)  and  0.1°  (bottom) 

Comparison  of  the  resulting  lift  and  drag  coefficients  (CL  and  CD,  respectively) 
obtained  from  the  quasi-steady  sweep  procedure  to  those  obtained  from  the  steady- 
state  methodology  were  reasonable  for  both  da  values.  For  the  pitching  moment, 
Cm,  the  results  obtained  from  the  sweep  procedure  showed  reasonably  good 
agreement  between  the  2  da  increments  investigated.  However,  there  are  some  very 
small  differences  between  the  steady-state  values  and  the  sweep  values  above  a  = 
5°  (Fig.  5)  that  are  most  noticeable  in  the  fin  component  contribution.  For  the  roll 
torque,  Cj  ,  there  was  good  agreement  between  the  steady-state  results  and  the 

sweep  procedure  for  da  =  0.1°;  agreement  is  not  as  good  for  the  sweep  procedure 
with  da  =  0.25°  (Fig.  6).  The  da  =  0.25°  sweep  results  are  not  symmetric  about 
a  =  0°  as  is  expected  for  this  symmetric  configuration.  Therefore,  the  convergence 
of  the  inner  iterations  for  da  =  0.25°  was  investigated. 
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Fig.  5  Pitching-moment  coefficient  for  BF  comparing  da  variation  to  steady-state  results, 
CFL  =  10  and  N  =  50 


0.098 

0.097 

0.096 

0.095 

0.094 

_o 

^  0.093 
0.092 
0.091 
0.090 
0.089 
0.088 


a 


Fig.  6  Roll-torque  coefficient  for  BF  comparing  da  variation  to  steady-state  results,  CFL 
=  10  and  N  =  50 
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As  expected,  increasing  the  number  of  inner  iterations  allowed  further  convergence 
in  the  solution  residuals  for  the  inner  iterations  (Fig.  7).  Nearly  2-orders-of- 
magnitude  drop  are  achieved  for  200  inner  iterations.  With  the  improved 
convergence  at  each  time  step,  a  change  in  Qo  over  the  first  40  steps  of  the  sweep 

procedure  can  be  seen  (Fig.  8).  Although  the  actual  values  do  not  change 
significantly  (the  scale  is  expanded),  the  shift  to  the  left  indicates  a  better  solution 
can  be  achieved  (similar  to  that  of  da  =  0.1°)  with  an  increase  in  the  number  of 
sub  iterations. 


Residual  Plot  (Normalized) 


Global  Iterations 


Residual  Plot  (Normalized) 


Global  iterations 


Fig.  7  Time  history  of  solution  residual  as  a  function  of  number  of  inner  iterations  for  CFL 

=  10,  da  =  0.25°;  N  =  100  (top)  and  N  =  200  (bottom) 
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Fig.  8  Roll-torque  coefficient  for  BF  for  varying  number  of  inner  iterations,  CFL  =  10  and 
da  =  0.25° 


Using  da  =  0.1°  and  N  =  50  or  da  =  0.25°  and  N  =  200  both  produce  reasonable 
sweep  results  when  compared  with  the  steady-state  solutions,  in  a  reasonably 
efficient  manner.  All  solutions  were  obtained  using  256  central  processing  units 
(CPUs)  on  a  SGI  ICE  X  supercomputer.  Each  steady-state  solution  takes 
approximately  3.5  wall-clock  h  (2,000  iterations  to  converge,  approximately  6  s  per 
iteration).  For  da  =  0.25°  and  N  =  200,  it  takes  approximately  2 1  min  per  angle  or 
40  h  for  the  entire  sweep  (10,000  CPU  h).  For  da  =  0.1°  and  N  =  50,  this  time  was 
further  reduced  to  3.5  min  per  time-step  or  15  h  for  the  entire  sweep  (4,000  CPU 
h).  Even  assuming  only  1°  increments,  it  would  take  nearly  26,000  CPU  h  for  the 
29  steady-state  solutions,  thus  realizing  a  CPU  savings  of  up  to  85%  using  the 
quasi-steady  sweeps. 

4.1.2  CFL  =  50 

After  completing  the  sweep-procedure  comparison  for  the  BF  with  CFL  =  10,  it 
was  determined  the  quasi-steady  sweep  procedure  could  be  completed  in  a  more 
efficient  manner  if  a  larger  CFL  value — similar  to  that  used  for  the  steady-state 
simulations — was  chosen,  rather  than  using  a  small  CFL  as  would  be  the  case  if  a 
truly  time-accurate  solution  was  being  sought.  Use  of  a  larger  CFL  in  a  steady-state 
simulation  typically  allow  convergence  in  fewer  iterations.  Since  a  quasi-steady 
sweep  is  solving  for  a  “steady-state”  solution  at  each  time-step,  it  was  thought  that 
convergence  at  each  time-step  could  be  achieved  with  fewer  inner  iterations  if  a 
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larger  value  for  the  CFL  number  was  chosen.  As  the  CFL  number  for  some  of  the 
steady-state  simulations  at  higher  angles  of  attack  had  to  be  limited  to  50  in  order 
for  convergence  to  be  achieved  for  the  BF  (if  a  CFL  greater  than  50  was  chosen, 
convergence  was  not  possible),  it  was  desirable  to  know  the  improved  efficiency  of 
sweep  procedure  and  the  required  parameters  if  the  CFL  was  set  to  50  for  the  quasi¬ 
steady  sweeps. 

Because  inner-iteration  convergence  is  required  to  achieve  accurate  solutions,  the 
number  of  inner  iterations  required  for  convergence  at  each  da  (0.5°,  0.25°,  and 
0.1°)  was  investigated  first.  The  findings  for  each  da  were  the  same;  there  were 
small  differences  in  force  and  moment  coefficients  when  the  results  of  N  =  50  and 
N  =  100,  but  virtually  no  differences  when  the  results  of  N  =100  and  N  =  200  were 
compared. 

The  convergence  history  of  the  inner  iterations  for  the  da  =  0.25°  case  (Fig.  9) 
shows  why  this  is  the  case.  There  is  a  decrease  in  the  residuals  during  the  inner 
iterations  from  approximately  2E-2  to  IE-2  when  N  =  100  is  used  rather  than 
N  =  50.  When  N  is  further  increased  to  200,  the  inner-iteration  residuals  plateau  at 
IE-2  indicating  that  the  solution  has  converged.  Thus,  N  =  100  was  chosen  for  all 
da  investigated. 
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Residual  Plot  (Normalized) 
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Residual  Plot  (Normalized) 


■energy(ESlOl) 
mass(ES_1  01 ) 
■x-momentum(ES_101) 
y-momentum(ES_1  01) 
■  z-momentum(ES_1  01 ) 


Global  iterations 


Fig.  9  Time  history  of  solution  residual  as  a  function  of  number  of  inner  iterations  for  CFL 
=  50,  da  =  0.25°;  N  =  50  (top),  N  =  100  (middle),  and  N  =  200  (bottom) 


Approved  for  public  release;  distribution  is  unlimited. 


15 


With  inner  iteration  convergence  achieved,  the  results  of  the  quasi-steady  sweep 
for  da  =  0.5°,  da  =  0.25°  and  da  =  0.1°  were  compared.  Most  of  the  forces  and 
moments  for  the  total  configuration,  as  well  as  the  body  component  contribution 
and  fin  component  contribution  (C^)  were  identical  and  in  good  agreement  with  the 
steady-state  solutions.  Figures  10-13  show  the  fin  component  contribution  to  the 
forces  and  moments  for  the  sweep  procedure.  Only  Cio  shows  any  difference  for 

da  =  0.5°,  and  that  is  negligible  even  with  the  enlarged  scale  used.  The  total  value 
of  roll  torque  shown  in  Fig.  13  is  equivalent  to  the  fin  component  contribution  as  it 
is  the  only  contribution. 
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Fig.  10  Fin  component  contribution  to  drag  coefficient  for  BF  with  CFL  =  50,  N  =  100,  da 
comparison 
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Fig.  11  Fin  component  contribution  to  lift  coefficient  for  BF  with  CFL  =  50,  N  =  100,  da 
comparison 


Fig.  12  Roll-torque  coefficient  for  BF  with  CFL  =  50,  N  =  100,  da  comparison 
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Fig.  13  Fin  component  contribution  to  pitching-moment  coefficient  for  BF  configuration 
with  CFL  =  50,  N  =  100,  da  comparison 

Figures  14-16  show  the  body  component  contribution  (CB)  to  the  forces  and 
moments  for  the  sweep  procedure.  Roll-torque  coefficient  is  not  shown  as  the  body 
does  not  contribute  to  that  moment.  Even  using  da=0.5°,  very  good  agreement  with 
the  steady-state  solutions  is  achieved.  Thus,  agreement  in  the  total  forces  and 
moments  predicted  using  the  sweep  procedure  and  the  steady-state  solutions  is 
achieved  for  all  of  the  increments  investigated,  indicating  that  for  the  BF  at 
subsonic  velocity  an  increment  as  large  as  da=0.5°  is  acceptable.  This  does,  of 
course,  assume  a  sufficient  number  of  inner  iterations  is  utilized.  It  becomes  the 
researcher’s  choice  of  the  desired  increment,  based  on  efficiency  or  resolution. 
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Fig.  14  Body  component  contribution  to  drag  coefficient  for  BF  with  CFL  =  50,  N  =  100, 
da  comparison 


Fig.  15  Body  component  contribution  to  lift  coefficient  for  BF  with  CFL  =  50,  N  =  100,  da 
comparison 
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Fig.  16  Body  component  contribution  to  pitching-moment  coefficient  for  BF  with  CFL 
=  50,  N  =  100,  da  comparison 

For  da  =  0.5°,  N  =  100,  1 1.7  wall-clock  h  on  256  cores  of  the  SGI  ICE  X  were 
required  (2,995.2  CPU  h  total).  For  da  =  0.25°,  N=  100,  21.7  wall-clock  h  on  256 
cores  of  the  SGI  ICE  X  were  required  (5,555.2  CPU  h  total).  For  da  =  0.1°,  N  = 
1 00,  5 1 .9  wall-clock  h  on  256  cores  of  the  SGI  ICE  X  were  required  (1 3,286.4  CPU 
h  total). 

For  the  body-fin  configuration,  completing  the  quasi-steady  sweep  allows  a  range 
of  angles  of  attack  to  be  investigated  in  significantly  less  time  than  if  individual 
steady-state  simulations  were  used  (see  Table  1).  When  completing  the  quasi¬ 
steady  sweep  simulations,  it  is  critical  to  ensure  the  inner  iterations  are  converged. 
The  number  of  inner  iterations  required  to  achieve  convergence  varies  with  CFL 
number  as  well  as  Mach  number.  The  study  conducted  by  Sahu  and  Heavy  (2010) 
found  only  25-50  inner  iterations  were  required  to  achieve  time-step  convergence 
in  the  transonic  and  supersonic  regimes.  In  the  subsonic  regime  investigated  here, 
50-200  inner  iterations  were  required  depending  on  the  CFL  number  chosen.  For 
CFL  =10,  the  required  number  of  inner  iterations  for  convergence  varied  with  da. 
Only  50  inner  iterations  were  required  for  da  =  0.10°,  while  200  inner  iterations 
were  required  for  da  =  0.25°.  This  meant  the  coarser  angular  resolution  actually 
required  more  computational  resources  (Table.  Additionally,  inner-iteration 
convergence  would  need  to  be  investigated  each  time  a  different  angular  resolution 
was  desired  if  CFL  =10  was  chosen.  For  CFL  =  50,  the  required  number  of  inner 
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iterations  for  convergence  (N  =  100)  did  not  vary  with  da.  Choosing  CFL  =  50 
ensures  that  once  inner-iteration  convergence  is  achieved  at  one  angular  resolution, 
changing  the  angular  resolution  does  not  vary  those  criteria.  The  computational 
requirement  is  directly  proportional  to  your  desired  angular  resolution. 

Table  1  Comparison  of  time  required  to  obtain  converged  solution  for  -14°  <  a  <  14°  at  da 
indicated 


CFL 

da  (°) 

N 

CPU  hours 

Steady-state 

1.00 

26,000 

10 

0.10 

50 

4,000 

0.25 

200 

10,000 

50 

0.10 

100 

13,286 

0.25 

100 

5,555 

0.50 

100 

2,995 

4.2  Body-Fin-Canard 

The  body-fin-canard  configuration  is  investigated  for  canards  in  the  x-orientation 
only  (i.e.,  45°  rotation  out  of  the  lift  plane).  The  canards  are  considered  in  their 
undeflected  position  ( 8  =  0°)  as  well  as  for  a  4°  deflection  about  its  quarter  chord 
for  a  pitch  up  maneuver  ( 8  =  4°).  Two  canard  deflection  angles  were  considered 
as  it  was  unknown  if  variations  in  canard  separation,  canard  stall,  canard-root 
and/or  -tip  vortex  shedding,  and  vortex-fin  interactions  due  to  canard  deflection 
would  require  different  values  of  the  critical  parameters  of  the  sweep  procedure  to 
obtain  agreement  with  the  steady-state  solutions. 

As  was  done  for  the  BF,  steady-state  simulations  were  completed  first.  The  steady- 
state  simulations  were  completed  at  a  =  0°,  +2°,  ±3°,  ±5°,  +8°,  +10°,  and  + 
14°.  The  sweep  procedure  was  performed  using  the  steady-state  solution  at  a  = 
—  14°  as  the  initial  condition.  The  critical  parameters  were  varied  to  determine  the 
requirements  for  an  accurate  solution.  The  critical  parameters  varied  were  CFL  = 
10  and  50;  da  =  0.25°,  0.1°,  and  0.05°;  andN  =  50,  100,  and  200. 

4.2.1  S=  0° 

As  with  the  BF,  the  body-tail-canard  configuration  with  the  canards  undeflected 
(BFC0),  CFL  =  10,  N  =  50  was  investigated  first  for  da  =  0.25°  and  da  =  0.1°. 
CFL  =10  was  considered  for  the  BFC0  as  the  BFC0  quasi-steady  sweeps  were 
completed  simultaneously  with  the  quasi-steady  sweeps  for  BF.  The  CFL  =10 
results  are  presented  for  the  BFC0  to  further  show  the  effects  of  not  converging  the 
inner  iterations. 

Inner-iteration  convergence  for  the  BFC0  was  similar  to  that  of  the  BF  (Fig.  4). 
Agreement  between  the  quasi-steady  sweep  procedure  and  steady-state  results  for 
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the  BFCO  for  CFL  =10  were  not  as  good  as  those  of  the  BF,  indicating  inner- 
iteration  convergence  is  even  more  important  as  the  flow  field  becomes  more 
complicated.  Agreement  was  good  for  the  lift  and  drag  coefficients  (Figs.  17  and 
18,  respectively).  Flowever,  this  is  not  the  case  for  either  Cm  (Fig.  19)  or  Qo  (Fig. 

20);  agreement  was  quite  poor  for  |a|  >  5°,  especially  for  da  =  0.25°.  While  the 
discrepancies  in  Qo  were  also  present  in  the  BF,  those  in  Cm  were  not.  An 

examination  of  the  component  contributions  to  Cm  shows  the  canard  component 
contribution  ( Cc )  is  the  cause  of  the  discrepancy  (Fig.  21).  With  the  inner  iterations 
not  adequately  converged,  canard  stall  and  separation  are  not  properly  predicted. 
Increasing  N  =  200  for  da  =  0.25°  produces  results  closer  to  those  for  da  =  0.1°, 
N  =  50  (Fig.  22),  but  at  a  significant  time  penalty. 


Fig.  17  Lift  coefficient  for  BFCO  with  CFL  =  10,  N  =  50,  da  comparison 
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Fig.  18  Drag  coefficient  for  BFCO  with  CFL  =  10,  N  =  50,  da  comparison 


Fig.  19  Pitching-moment  coefficient  for  BFCO  with  CFL  =  10,  N  =  50,  da  comparison 
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Fig.  20  Roll-torque  coefficient  for  BFCO  with  CFL  =  10,  N  =  50,  da  comparison 


Fig.  21  Component  contributions  to  pitching-moment  coefficient  for  BFCO  with  CFL  =  10, 
N  =  50,  da  comparison 


Approved  for  public  release;  distribution  is  unlimited. 

24 


Fig.  22  Improved  agreement  for  pitching-moment  and  roll-torque  coefficients  for  BFCO 
for  CFL  =  10,  with  N  =  200  for  larger  da 

While  it  is  possible  that  a  smaller  da  or  an  increased  number  of  inner  iterations 
could  lead  to  a  more  accurate  sweep  solution  for  the  BFCO,  the  amount  of 
computing  time  required  to  achieve  the  more  accurate  solutions  would  eliminate 
the  efficiency  of  the  quasi-steady  sweep  procedure.  Thus,  da  =  0.1  °,N  =  100  was 
not  investigated  for  improved  agreement  between  the  static  and  quasi-steady  sweep 
results  for  Cm  and  Cio.  Based  on  the  results  of  the  BF — using  a  larger  CFL,  similar 

to  that  used  to  converge  the  steady-state  solutions,  for  the  quasi-steady  sweeps 
produced  more  reliable  results — the  remainder  of  the  quasi-steady  sweeps  for  the 
BFCO  were  completed  with  CFL  =  50. 

As  was  the  case  for  quasi-steady  sweeps  on  the  BF  when  CFL  =  50,  N  =  100  was 
found  to  produce  a  converged  solution  for  both  da  investigated  (0.25°,  0.1°),  the 
size  of  da  did  not  matter  with  the  larger  CFL.  The  inner- iteration  convergence 
histories  were  similar  to  those  found  for  the  BF  (Fig.  9).  For  CFL  =  50,  N  =  100, 
significantly  better  agreement  was  found  between  the  static  results  and  the  quasi¬ 
steady  sweep  results  for  both  Cm  (Fig.  23)  and  Cto  (Fig.  24)  for  the  BFCO.  Prediction 

of  the  body  and  fin  component  contributions  to  Cm  by  the  sweep  procedure 
remained  good  with  the  better-converged  inner  iterations  as  would  be  expected. 
Canard  stall  was  now  properly  captured  (when  compared  to  the  steady-state 
solutions)  with  proper  inner-iteration  convergence  even  for  da  =  0.25°  (Fig.  25). 
While  it  is  possible  an  adequately  converged  sweep  solution  could  still  be  obtained 
for  da  =  0.5°,  it  was  decided  that  da  =  0.25°  was  sufficiently  large  enough  due 
to  the  small  discrepancies  that  were  being  observed  in  the  BF. 
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Fig.  24  Roll-torque  coefficient  for  BFCO  with  CFL  =  50,  N  =  100,  da  comparison 
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Fig.  25  Single  canard  component  contribution  to  pitching-moment  coefficient  for  BFCO 
with  CFL  =  50,  N  =  100,  da  comparison 

The  time  savings  achieved  in  obtaining  the  solutions  for  the  BFCO  using  the  quasi¬ 
steady  sweep  procedure  as  compared  to  the  steady-state  method  is  even  greater  than 
it  was  for  the  BF.  The  BFCO  solutions  were  completed  on  384  cores  of  the  SGI  ICE 
X  supercomputer  SPIRIT  at  the  US  Air  Force  Research  Laboratory’s  Department 
of  Defense  Supercomputing  Resource  Center.  A  single  steady-state  solution  could 
be  obtained  in  approximately  5  h.  That  means,  just  to  obtain  the  1 3  steady-state 
solutions  for  comparison  here  took  65  wall-clock  h  or  24,960  CPU  h.  The  sweep 
procedure  using  da  =  0.25°  (CFL  =  50,  N  =  100)  took  just  37.7  wall-clock  h 
(14,477  CPU  h),  a  42%  time  savings.  This  would  significantly  improve  resolution 
for  predicting  stall  while  greatly  reducing  resources.  Increasing  resolution  further 
by  using  da  =  0.1°  (CFL  =  50,  N  =  100)  increases  use  of  computing  resources  to 
37,900  CPU  h.  While  this  is  50%  increase  in  use  of  resources  over  that  used  to 
generate  the  steady-state  solutions  presented  here,  a  lot  more  data  were  generated 
with  those  resources. 


4.2.2  S=  4° 

The  final  case  considered  as  part  of  this  study  was  the  body-fin-canard 
configuration  with  the  canards  deflected  to  4°  for  pitch  control  (BFC4).  A 
deflected-canard  configuration  was  investigated  to  ensure  the  variations  in  the  flow 
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that  occur  when  the  canards  are  deflected  (separation,  stall,  vortex  shedding, 
vortex-fin  interaction,  etc.)  did  not  cause  a  change  in  sweep’s  critical  parameters. 

Convergence  criteria  for  CFL  =  1 0  were  again  investigated  because  it  was  thought 
that  a  more  complicated  configuration  (i.e.,  more  flow  nonlinearities)  may  require 
a  smaller  CFL  for  convergence  of  the  steady-state  solution  at  high  a ,  which  would 
in  turn  limit  the  CFL  for  the  sweep  procedure.  Due  to  the  simpler  BF  requiring  N 
=  200  for  convergence  of  the  quasi-steady  sweep  when  da  =  0.25°,  only 
da  =  0.25°  withN  =  200  is  compared  to  the  results  for  da  =  0.1°  and  da  =  0.05° 
each  with  N  =  50,  which  was  the  required  number  of  inner  iterations  for 
convergence  on  the  BF. 

As  with  the  previous  configurations  considered,  only  Cm  and  C1q  showed  any 
noticeable  variation  with  da;  the  forces  converge  more  quickly  than  the  moments. 
The  variation  in  Cm  was  only  noticeable  at  |  a  \  >  10°  and  even  there  it  was  typically 
less  than  5%  (Fig.  26).  The  variation  in  Cio  was  only  noticeable  for 

\a\  >  8°  and  on  an  enlarged  scale  (Fig.  27).  This  further  validates  that  it  is  not  the 
da  that  matters  but,  rather,  that  a  sufficient  number  of  inner  iterations  are  completed 
such  that  convergence  is  reached  (Fig.  27)  for  the  chosen  CFL  number.  A  smaller 
da  requires  fewer  inner  iterations  and  shows  less  inner-iteration  residual  drop  as 
there  are  fewer  changes  in  the  flow  between  steps.  This  implies  that  if  there  are 
regions  of  a  where  the  flow  is  more  difficult  to  converge — thus  requiring  a  smaller 
CFL  number — the  quasi-steady-sweep  procedure  can  be  completed  with  an 
increased  number  of  inner  iterations  in  the  region;  a  smaller  da  may  not  be 
required.  This  would,  of  course,  require  additional  CPU  time  for  the  increased  inner 
iterations. 
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Fig.  26  Cm  comparison  of  BFC4  for  CFL  =  10  and  da  =  0.25°  with  N  =  200  to  da  =  0.1°  and 
da  =  0.05°  with  N  =  50 


a  (deg) 


Fig.  27  Cio  comparison  of  BFC4  for  CFL  =  10  and  da  =  0.25°  with  N  =  200  to  da  =  0.1°  and 
da  =  0.05°  with  N  =  50 
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Fig.  28  Inner-iteration  convergence  history  for  BFC4  for  CFL  =  10,  da=  0.25°  with  N  =  200 
(top),  da  =  0.1°  with  N  =  50  (middle),  and  da  =  0.05°  with  N  =  50  (middle) 
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Since  CFL  =  50,  N  =  100  produced  better  results  at  larger  da  than  a  CFL  =  10  for 
the  BFCO,  it  was  investigated  for  the  BFC4  to  determine  if  the  quasi-steady  sweep 
could  be  obtained  any  more  efficiently.  For  CFL  =  50,  convergence  of  the  inner- 
iteration  residuals  was  confirmed  for  N  =  1 00.  Residual  drops  of  approximately  2E- 
2  and  4E-2  were  achieved  for  da  =  0.25°  and  da  =  0.1°,  respectively  (Fig.  29). 
The  residual  drop  is  not  as  large  for  da  =  0.1°  because  there  is  less  change  in  the 
flow  due  to  the  smaller  change  in  a.  For  both  increments,  the  residuals  appear  to 
be  reaching  a  plateau,  which  indicates  convergence.  Therefore,  N  =  100  was  used 
and  no  further  inner-iteration  convergence  was  investigated  for  CFL  =  50. 


Residual  Plot  (Normalized) 


Global  iterations 


Residual  Plot  (Normalized) 


Global  iterations 


Fig.  29  Inner-iteration  convergence  history  for  BFC4  with  CFL  =  50,  N  =  100;  da  =  0.25° 
(top)  and  da=  0.1°  (bottom) 
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The  true  test  of  inner-iteration  convergence  is  the  comparison  of  the  quasi-steady- 
sweep  procedure’s  results  to  the  results  obtained  from  the  steady-state 
methodology.  All  force  and  moment  coefficients  (Figs.  30-33)  show  very  good 
agreement  with  the  solutions  from  the  steady-state  methodology,  indicating  the 
inner  iterations  must  indeed  be  converged.  Even  the  nonlinearity  in  the  moment 
coefficients  is  accurately  captured.  From  the  analysis  of  C £  for  each  canard 
(Fig.  34),  canard-stall  angle  is  known  to  within  a  0.25°;  no  interpolation  is 
necessary. 


Fig.  30  Drag  coefficient  for  BFC4  with  CFL  =  50,  N  =  100,  da  comparison 
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Fig.  31  Lift  coefficient  for  BFC4  with  CFL  =  50,  N  =  100,  da  comparison 


Fig.  32  Roll-torque  coefficient  for  BFC4  with  CFL  =  50,  N=100,  da  comparison 
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Fig.  34.  Individual  canards’  pitching-moment  contributions  with  da  =  0.25°,  CFL  =  50, 
N  =  100 
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As  the  mesh  size  for  the  BFC4  is  approximately  the  same  as  that  of  the  BFCO,  the 
required  computational  time  was  found  to  be  similar.  For  da  =  0.25°,  36.5  wall- 
clock  h  on  384  cores  of  the  SGI  ICE  X  cluster  (approximately  14,000  CPU  h)  were 
required.  For  da  =  0.1°,  86.3  wall-clock  h  on  384  cores  of  the  SGI  ICE  X  cluster 
(approximately  33,000  CPU  h)  were  necessary.  Compare  this  to  a  required  wall- 
clock  time  of  3.8  h  (1,460  CPU  h)  for  the  steady-state  solution.  The  break-even 
point  is  10  steady-state  solutions.  If  more  than  10  steady-state  solutions  are  required 
to  adequately  resolve  the  static  coefficients,  which  is  likely  the  case  of  a  highly 
nonlinear  flow  such  as  this,  the  sweep  procedure  is  definitely  the  more  efficient 
method.  The  increment  is  dependent  on  the  resolution  desired  as  was  the  case  with 
the  BFCO.  The  physics  is  being  adequately  captured. 

5.  Conclusions 


A  quasi-steady  sweep  procedure  that  allows  determination  of  the  static 
aerodynamic  coefficients  for  a  range  of  angles  of  attack  in  a  single  simulation  is 
extended  to  subsonic  Mach  numbers.  Given  a  steady-state  solution,  the  procedure 
can  be  used  to  generate  multiple  quasi-steady  numerical  simulations  over  a  select 
range  of  angles  of  attack  for  a  given  set  of  flight  conditions.  Computational  fluid- 
dynamics  simulations  are  completed  for  the  ARL-designed  HMA,  a  canard- 
controlled,  fin-stabilized  projectile,  at  Mach  0.65  for  a  moderate  range  of  angles  of 
attack  (—14°  <  a  <  14°)  using  the  sweep  procedure.  Separate  steady-state 
solutions  were  also  computed  at  various  angles  of  attack  within  this  range  for 
comparison. 

Three  configurations  of  the  HMA  geometry  were  considered  to  evaluate  the  sweep 
procedure  for  increasingly  complex  geometries:  1)  body-fin,  2)  body-fin-canard 
at  0°  canard-deflection  angle,  and  3)  body-fin-canard  with  canards  deflected  4°  for 
a  pitch  up  maneuver.  For  each  configuration,  the  parameters  of  da,  CFL  number, 
and  number  of  inner  iterations  were  varied  to  find  the  most  efficient  combination 
to  achieve  a  converged  solution.  The  initial  values  used  for  these  parameters,  based 
on  previous  findings  for  supersonic  velocities,  were  found  to  be  insufficient  in  the 
subsonic  regime  to  find  agreement  with  the  steady-state  solution  for  Cm  and  Cio. 

Inner-iteration  convergence  is  critical  to  accurate  solutions  for  the  quasi-steady 
sweep  procedure.  Depending  on  the  CFL  number  and  da  chosen,  inner-iteration 
convergence  is  achieved  when  a  magnitude  drop  of  greater  than  IE-2  in  residuals 
occurs  or  the  residual  levels  do  not  vary  significantly;  the  inner-iteration 
convergence  needs  to  be  checked  periodically  during  the  sweep  to  ensure  there  is 
not  an  angle  of  attack  that  is  having  difficulty  converging.  If  inner  iterations  are 
having  difficulty  converging  during  just  a  section  of  the  sweep  trajectory,  the  CFL 


Approved  for  public  release;  distribution  is  unlimited. 

35 


number  could  be  reduced  and  either  the  da  could  be  reduced  or  the  number  of  inner 
iterations  increased  (or  both)  for  that  section  to  ensure  convergence  throughout  the 
range  of  angles  of  attack. 

Variation  of  CFL  number  in  the  sweep  procedure  showed  that  agreement  with  the 
steady-state  solution  could  be  approached  for  either  CFL  =  10  or  CFL  =  50  for  any 
of  the  3  configurations  if  a  sufficient  number  of  inner  iterations  and  a  sufficiently 
small  da  were  used.  Using  CFL  =  50,  da  =  0.25°,  and  N  =  100,  the  sweep 
procedure  was  found  to  achieve  excellent  agreement  with  the  aerodynamic  forces 
and  moments  obtained  from  the  steady-state  method  for  all  3  configurations.  The 
CFL  =10  required  the  use  of  a  smaller  da  (0.1°  versus  0.25°)  and/or  an  increased 
number  of  inner  iterations  (200  versus  100)  to  achieve  convergence  as  compared  to 
CFL  =  50  to  obtain  results  that  were  in  reasonable  agreement  with  the  Cm  and  Q 

values  obtained  from  the  steady-state  method.  Thus,  convergence  occurs  more 
quickly  when  a  larger  CFL  number  is  used,  as  would  be  expected.  However,  if  the 
CFL  number  needs  to  be  limited  to  achieve  a  converged  flow  field,  inner-iteration 
convergence  can  still  be  achieved  with  use  of  smaller  da  or  an  increased  number 
of  inner  iterations. 

Although  smaller  da  values  and  more  inner  iterations  are  required  to  obtain 
accurate  moment  coefficients  from  the  sweep  procedure  for  a  subsonic  flow,  the 
quasi-steady  sweep  procedure  remains  an  efficient  means  of  obtaining  a  range  of 
angles  of  attack  for  a  given  set  of  flow  conditions.  A  typical  steady-state  solution 
at  one  angle  of  attack  for  the  BFC4  took  1 ,460  CPU  h  (approximately  3.5  h  on  384 
cores).  The  sweep  procedure  with  CFL  =  50,  da  =  0.25°,  N  =  100  required 
approximately  14,000  CPU  h  (approximately  36.5  h  on  384  cores).  This  means  that 
in  the  time  it  would  take  to  complete  10  steady-state  solutions,  a  solution  with  112 
angles  was  generated.  For  a  complex  configuration  that  may  potentially  have 
nonlinear  forces  and  moments,  this  increased  resolution  is  critical.  For  a  less 
complex  configuration  such  as  the  BF,  where  the  forces  and  moments  are  known 
to  be  linear,  the  increased  resolution  of  the  sweep  procedure  is  not  necessary; 
greater  efficiency  could  be  achieved  by  just  choosing  a  small  number  of  angles  of 
attack  to  simulate  with  the  steady-state  method  and  using  standard  curve  fit 
procedures  to  generate  data  at  the  other  angles. 
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Appendix.  Quasi-Steady  Sweep's  Implementation  in  CFD++ 
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The  CFD++  implementation  of  the  sweep  procedure  is  fairly  typical  of  any 
computational  fluid  dynamics  solver.  There  are,  however,  a  number  of  nuances  in 
this  particular  solver  that  the  author  believed  should  be  documented.  All  have  to  do 
with  the  setup  of  the  procedure  and  not  the  actual  solution. 

In  initializing  the  quasi-steady  sweep  procedure  (and  all  time-accurate  simulations) 
in  CFD 1 1 ,  one  must  first  determine  how  they  want  to  start:  reset  the  step  to  zero  or 
continue  from  the  last  step  in  the  steady-state  solution;  the  time  always  starts  from 
a  value  of  0.  If  one  chooses  to  continue  from  the  last  step  in  the  steady-state 
solution,  the  following  does  not  affect  the  solution;  but,  if  one  chooses  to  reset  the 
step  to  zero — the  author’s  preference  for  the  sweep  procedure — the  following  must 
be  completed.  There  are  a  number  of  equations  that  are  not  turned  on  at  Step  0  by 
default.  These  include  preconditioning  and  turbulence  modeling.  These  can  be 
easily  set  to  start  at  Step  0  either  in  the  steady  to  unsteady  wizard  or  in  the  Reimann 
Solver  and  Turbulence  Control  panels,  respectively.  The  author  recently  found  the 
additional  transport  equations  are  solved  in  nonconservation  mode  for  the  first  20 
steps  in  CFD 1 1  by  default.  Prior  to  Version  15.1  Update  4,  one  had  to  add  the  option 
“iftcon_ntstart  0”  to  the  mcfifinp  input  file  such  that  the  additional  transport 
equations  are  solved  in  conservation  mode  from  Step  0.  In  Version  15.1  Update  4 
and  later,  there  is  an  “Extra  Eqns”  button  in  the  second  panel  of  the  Time  Integration 
Panel.  By  choosing  this  “Extra  Eqns”  Panel,  one  can  change  the  conservation 
updates  for  the  equations  for  turbulence,  scalars,  and  species  to  start  at  Step  0.  This 
ensures  that  every  single  step  is  physically  relevant.  Although  not  relevant  to  this 
subsonic  simulation,  one  also  needs  to  ensure  that  1st-  to  2nd-order  spatial 
discretization  blending  is  turned  off  and  that  the  Right  Hand  Side  dissipation  is  0 
for  the  time-accurate  solution  that  starts  at  Step  0  if  either  was  used  to  obtain  the 
steady-state  simulation. 

CFD++  uses  file-based  rotation  to  complete  the  mesh  motion  used  for  the  sweep 
procedure.  For  the  case  considered  here,  only  single-axis  rotational  motion  in  the 
body  frame  is  required.  As  per  the  CFD++  help,  the  file  format  is  as  follows: 

Origin  (x  y  z  values) 

Direction  vector  (x  y  z  values) 

#  datalines 

Timel  (s)  rotation  anglel  (radians)  rotation  ratel  (radians/s) 

Time2  (s)  rotation_angle2  (radians)  rotation_rate2  (radians/s) 
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For  an  angle-of-attack  sweep,  the  origin  should  be  at  the  center  of  gravity.  If  the 
origin  of  the  computational  domain  is  at  the  center  of  gravity,  this  is  (0,  0,  0).  If  the 
origin  of  the  computational  domain  is  at  another  location,  the  center  of  gravity’s 
location  relative  to  the  origin  of  the  computational  domain  should  be  entered. 

The  directional  vector  is  self-explanatory.  For  an  alpha  sweep,  with  the  x-axis 
aligned  with  (and  pointing  downstream)  the  flow  and  the  z-axis  up,  the  rotational 
axis  is  the  y-axis.  Therefore  the  directional  vector  would  be  (0  1  0). 

For  an  angle-of-attack  sweep  where  a  single  angle-of-attack  increment  (da)  is 
utilized,  only  2  data  lines  are  required.  Because  the  sweep  procedure  always  starts 
from  time  zero,  regardless  of  whether  or  not  you  reset  the  step  to  zero,  the  first  data 
line  is  at  time  zero  with  a  rotation  angle  of  zero.  The  rotation  angle  is  zero  because 
the  steady-state  solution  from  which  the  quasi-steady  sweep  is  initialized  is 
obtained  on  the  computational  domain  that  has  been  rotated  to  the  minimum  angle 
of  attack  desired.  The  second  data  line  is  a  bit  more  complicated.  The  final  time 
(Time finai)  is  the  total  angle  of  sweep  (asweep)  divided  by  the  increment  (Eq.  A- 
1).  The  rotation  angle  is  the  total  angle  of  sweep  in  radians,  which  for  —14°  <  a  < 
14°  is  28°  or  0.4886922  radians.  CFD++  calculates  the  time-accurate  solution  at  the 
starting  angle  of  attack  on  Step  1 ;  the  mesh  is  not  rotated  until  the  end  of  Step  1 . 
Therefore,  to  obtain  the  full  range  of  the  sweep  it  is  necessary  to  calculate  a  solution 
for  one  additional  increment.  This  means  a  total  of  1 1 3  (including  0)  seconds  (steps) 
are  simulated. 

Timefinal  =  ^  (A-l) 

The  contents  of  the  file  for  the  alpha  sweep  performed  here  is  shown  below: 

0  0  0 

0  1  0 

2 

0  0.0000000  0 

112  0.4886922  0 


If  more  than  one  da  needs  to  be  utilized  to  ensure  resolution,  an  additional  data  line 
can  be  added  for  each  additional  increment  desired.  The  first  data  line  remains  the 
same  as  in  the  previous  case  as  the  mesh  is  at  a  =  —14°  to  start  the  simulations. 
Equation  A-2  is  now  used  for  each  additional  data  line. 

Timei  =  +  Timei_1  ,  (A-2) 

where 
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a  sweep  i  ai  final  ^initial  •  (A-3) 

aifinai  's  ^1C  final  angle  in  the  increment  and  cciinitlal  is  the  initial  angle  on  the 
increment.  The  value  for  the  second  column  is  determined  by 

rotation_angle*  =  aifinal  -  ainitial  (A-4) 


where  0Cintiai  is  the  initial  angle  of  attack.  For  example,  if  it  was  desired  to  use 
da  =  0.1  for  —14°  <  a  <  —9°,  da  =  0.25  for  —9°  <  a  <  3°,  and  da  =  0.1  for 
3°  <  a  <  14°,  the  content  of  the  file  for  the  alpha-sweep  procedure  would  be  as 
shown  below  and  a  total  of 209  s  (i.e.,  steps)  simulated. 

0  0  0 

0  1  0 

4 

0  0.0000000  0 

50  0.0872665  0 

98  0.2967060  0 

208  0.4886922  0 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


3-D 

ARL 

BF 

BFCO 

BFC4 

Cd 

CL 

Ci 

Lo 

r 

da 

CFD 

CFL 

CPU 

HMA 

M 

N 

RANS 

RBD 

Rt 

y+ 

a 

k 

£ 

Superscripts 

B 

F 

C 


3 -dimensional 

US  Army  Research  Laboratory 

body-fin  configuration 

body-fin-canard  configuration  with  5  =  0° 

body-fin-canard  configuration  with  5  =  4° 

drag  coefficient 

lift  coefficient 

roll-torque  coefficient 

pitching-moment  coefficient 
angle-of-attack  increment 
Computational  Fluid  Dynamics 
Courant-Friedrichs-Lewy  number 
central  processing  unit 
high-maneuverability  airframe 
Mach  number 
number  of  inner  iterations 
Reynolds-Averaged  Navier-Stokes 
rigid-body  dynamics 
undamped  eddy  viscosity 
nondimensional  wall  distance 
angle  of  attack 
turbulent  kinetic  energy 
dissipation  rate 

body  component 
fin  component 
single-canard  component 
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